integration of 10x and inDrop data

library(Seurat)
library(Matrix)
library(tidyverse)
library(viridis)

Load the Seurat object

Before we load the object, we’re going to switck the rownames of our data to be the symbols. There is a 1-1 mapping of gene ID to symbol in these samples so we can do this with no loss of information.

dat = readRDS(file.path("data", "counts.rds"))
md = readRDS(file.path("data", "metadata.rds"))
g2s = readr::read_csv(file.path("metadata", "tx2gs.csv"),
                      col_names=c("gene", "symbol", "transcript")) %>%
  dplyr::select(-transcript) %>%
  unique()
rownames(dat) = g2s$symbol[match(rownames(dat), g2s$gene)]
seurat = CreateSeuratObject(dat, meta.data=md)

Overview

First we’ll look at the uncorrected data; we’ll just correct for mitochondrial percentage and nothing else. Then we’ll run PCA, clustering and look at the UMAP plots to get an idea of how much run-to-run and technology variation we are dealing with.

seurat = SCTransform(seurat, vars.to.regress = "pctMito", verbose = FALSE)
seurat = RunPCA(seurat, verbose = FALSE)
seurat = RunUMAP(seurat, dims = 1:30, verbose = FALSE)

p1 = DimPlot(seurat, group.by = "technology")
p1

A lot! We can see that the 10x and the inDrop data are mostly separate from each other, but there is some overlaps. We can see looking at the four samples that the four samples are very different from each other as well.

DimPlot(seurat, group.by = "sample")

So we are going to have some work to do trying to overlay these on top of each other. There are a few ways to do that. One is to just regress out the effect of the technology, the other is to do something more complex like CCA. We do need to do something about it though if we want to combine these two datasets. However, even without doing anything about the technological variation we can take a peek at the data and see if we can pick out some patterns.

uncorrected clustering

Here we’ll do a rough clustering and then look at the expression of our markers, without correcting for which technology or experiment they are from, to orient outselves.

clustering

ElbowPlot(seurat, ndims=50)

seurat = FindNeighbors(seurat, dims=1:25)
seurat = FindClusters(seurat, resolution=0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 10605
## Number of edges: 332518
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9401
## Number of communities: 25
## Elapsed time: 0 seconds

DimPlot(seurat)

cell type markers

markers = readr::read_csv("metadata/cell_type_markers.csv")
knitr::kable(markers)

cell_identity	genes
ISC	Dl
ISC	Smvt
ISC	CG10006
pISC	sna
ISCR2R5	polo
ISC/EB	esg
EB	E(spl)m3-HLH
EB	E(spl)malpha-BFM
EB	E(spl)mbeta-HLH
EB	klu
EE	pros
EE	7B2
EE_subtype	AstA
EE_subtype	AstC
EE_subtype	CCHa1
EE_subtype	CCHa2
EE_subtype	Tk
EE_subtype	NPF
EE_subtype	Orcokinin
EE_subtype	Dh31
aEC	nub
aEC	Myo31DF
mEC	nub
mEC	Myo31DF
pEC	nub
pEC	Myo31DF
aEC	betaTry
pEC	lambdaTry
mEC	lab
mEC	Vha100-4
cardia	Pgant4
copper	PGRP-LA
copper	PGRP-LC
copper	Apoltp
VM	Vn
LFC	Ilp3
LFC	PGRP-SC1a
LFC	PGRP-SC1b
LFC	PGRP-SD
iron	PGRP-SC2
iron	ZIP1(Zip42C.1)

Below we can see that cluster 3 and cluster 21 are likely ISC+EB or EB cells. Cluster 9, 17 and 20 are liekly some subtypes of EE cells. Groups 1, 4, 11, 14 and 23 might be aEC cells. Groups 5, 13 and 15 might be pEC cells. Group 8 looks like mECs. And group 15 looks like cardia cells.

Cluster 12 looks like it might be LFC cells. Cluster 7 might be iron or copper cells. Cluster 6 looks like it could be an EC cell, but we don’t know what compartment.

We’re just missing 2, 10, 18, 19, and 22. Not bad.

cluster.averages = AverageExpression(seurat, return.seurat=TRUE, add.ident="sample")
DoHeatmap(cluster.averages, features=head(markers, 20)$genes, size=2)

DoHeatmap(cluster.averages, features=tail(markers, 21)$genes, size=2)

clusters = seurat@meta.data$seurat_clusters
rough_celltype = case_when(
  clusters %in% c(3, 21) ~ "ISC/EB",
  clusters %in% c(9, 17, 20) ~ "EE",
  clusters %in% c(0, 1, 4, 11, 14, 23) ~ "aEC",
  clusters %in% c(5, 13, 16) ~ "pEC",
  clusters %in% c(8, 24) ~ "mEC",
  clusters %in% c(6) ~ "EC",
  clusters %in% c(15) ~ "cardia",
  clusters %in% c(12) ~ "LFC",
  clusters %in% c(7) ~ "iron/copper",
  TRUE ~ as.character(clusters))
seurat@meta.data$rough_celltype = rough_celltype
Idents(seurat) = seurat@meta.data$rough_celltype
roughcelltype.averages = AverageExpression(seurat, return.seurat=TRUE, add.ident="sample")
DoHeatmap(roughcelltype.averages, features=head(markers, 20)$genes, size=2)

DoHeatmap(roughcelltype.averages, features=tail(markers, 21)$genes, size=2)

p3 = DimPlot(seurat, group.by="rough_celltype", label=TRUE)
p3

It looks like we are able to identify the cells across both technologies, which is very good.

knitr::kable(seurat@meta.data %>%
             group_by(rough_celltype, technology) %>%
             summarise(n=n()))

rough_celltype	technology	n
10	10x	76
10	inDrop	323
18	10x	31
18	inDrop	214
19	inDrop	243
2	10x	100
2	inDrop	683
22	10x	10
22	inDrop	77
aEC	10x	913
aEC	inDrop	3142
cardia	10x	32
cardia	inDrop	225
EC	10x	186
EC	inDrop	316
EE	10x	306
EE	inDrop	551
iron/copper	10x	206
iron/copper	inDrop	266
ISC/EB	10x	141
ISC/EB	inDrop	643
LFC	10x	102
LFC	inDrop	221
mEC	10x	252
mEC	inDrop	226
pEC	10x	624
pEC	inDrop	496

unmarked clusters

This leaves us with clusters 2, 10, 18, 19 and 22 as unknown clusters.

seurat@meta.data %>%
  filter(is.na(rough_celltype)) %>%
  pull(seurat_clusters) %>%
  unique()

## factor(0)
## 25 Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ... 24

trypsin genes

tryp = c("deltaTry", "gammaTry", "alphaTry", "betaTry",
         "epsilonTry", "thetaTry", "kappaTry", "lambdaTry",
         "iotaTry", "etaTry", "zetaTry")
r1 = c("deltaTry", "gammaTry", "Muc68D")
r2 = c("deltaTry", "gammaTry", "alphaTry", "betaTry", "epsilonTry", "thetaTry")
r3 = c("thetaTry")
r4 = c("kappaTry", "lambdaTry", "iotaTry", "etaTry", "zetaTry")
r5 = c("iotaTry", "etaTry", "zetaTry")
crop = c("Cyp312a1", "Cyp4e3", "Spn27A", "spz", "pot")
compartments = unique(c(r1, r2, r3, r4, r5, crop))
FeaturePlot(seurat, features=r2)

FeaturePlot(seurat, features=r4)

FeaturePlot(seurat, features=r5)

FeaturePlot(seurat, features=crop)

DoHeatmap(cluster.averages, features=compartments, size=2)

It looks like group 12 and 15 might be from R1. Group 0, 1, 4, 11, 14, 19, 23 and 24 are from either R2 or R3. Group 7 looks like it might be from R3. Groups 5, 6 and 16 look like they are from R4. Group 18 looks like it is from the crop.

clusters = seurat@meta.data$seurat_clusters
compartment = case_when(
  clusters %in% c(12, 15) ~ "R1",
  clusters %in% c(0, 1, 4, 11, 14, 19, 23, 24) ~ "R2",
  clusters %in% c(7) ~ "R3",
  clusters %in% c(5, 6, 16) ~ "R4",
  clusters %in% c(18) ~ "crop",
  TRUE ~ as.character(clusters))
seurat@meta.data$compartment = compartment
Idents(seurat) = seurat@meta.data$compartment
compartment.averages = AverageExpression(seurat, return.seurat=TRUE, add.ident="sample")
DoHeatmap(compartment.averages, features=compartments, size=2)

p5 = DimPlot(seurat, group.by="compartment", label=TRUE)
p5

It looks like we are able to identify compartments cells are from across both technologies, which is very good.

knitr::kable(seurat@meta.data %>%
  group_by(compartment, technology) %>%
  summarise(n=n()))

compartment	technology	n
10	10x	76
10	inDrop	323
13	10x	161
13	inDrop	113
17	10x	165
17	inDrop	85
2	10x	100
2	inDrop	683
20	10x	35
20	inDrop	138
21	10x	140
21	inDrop	4
22	10x	10
22	inDrop	77
3	10x	1
3	inDrop	639
8	10x	252
8	inDrop	208
9	10x	106
9	inDrop	328
crop	10x	31
crop	inDrop	214
R1	10x	134
R1	inDrop	446
R2	10x	913
R2	inDrop	3403
R3	10x	206
R3	inDrop	266
R4	10x	649
R4	inDrop	699

write seurat

saveRDS(seurat, file.path("results", "seurat-rough.rds"))

integration

Here we wil use the integration functions in Seurat to integrate the 10x and the inDrop datasets.

Idents(seurat) = seurat@meta.data$seurat_clusters
indrop = SubsetData(seurat, subset.name="technology", accept.value="inDrop")
indrop = NormalizeData(object = indrop)
indrop = ScaleData(object = indrop)
indrop = FindVariableFeatures(indrop, do.plot = FALSE)
tenx = SubsetData(seurat, subset.name="technology", accept.value="10x")
tenx = NormalizeData(object = tenx)
tenx = ScaleData(object = tenx)
tenx = FindVariableFeatures(tenx, do.plot = FALSE)
hvf = union(VariableFeatures(indrop), VariableFeatures(tenx))
anchors = FindIntegrationAnchors(object.list=c(indrop, tenx), dims=1:30)
integrated = IntegrateData(anchorset = anchors, dims = 1:30, features.to.integrate=rownames(indrop))
DefaultAssay(integrated) = "integrated"

integrated clustering

Here we use the integrated data to re-do the clustering. After integration, we can see that the 10x and the inDrop datasets overlap with each other much better. We can also get a better idea about what the identities of cells might be.

integrated = ScaleData(integrated, verbose = FALSE)
integrated = RunPCA(integrated, npcs = 30, verbose = FALSE)
integrated = RunUMAP(integrated, reduction = "pca", dims = 1:30)
integrated = FindNeighbors(integrated, reduction = "pca", dims = 1:30)
integrated = FindClusters(integrated, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 10605
## Number of edges: 407093
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9284
## Number of communities: 22
## Elapsed time: 0 seconds

p2 = DimPlot(integrated, reduction = "umap", group.by = "technology")
p4 = DimPlot(integrated, reduction = "umap", group.by = "rough_celltype", label=TRUE)
p6 = DimPlot(integrated, group.by="compartment", label=TRUE)
library(gridExtra)
grid.arrange(p1 + ggtitle("before integration"), p2 + ggtitle("after integration"))

grid.arrange(p3 + ggtitle("before integration"), p4 + ggtitle("after integration"))

grid.arrange(p5 + ggtitle("before integration"), p6 + ggtitle("after integration"))

Here we redo the marker expression, using the integrated cluster IDs rather than the unintegrated IDs from before.

DimPlot(integrated, group.by="seurat_clusters", label=TRUE)

integrated.averages = AverageExpression(integrated, return.seurat=TRUE, add.ident="sample")
DoHeatmap(integrated.averages, features=head(markers$genes, 20), size=2)

DoHeatmap(integrated.averages, features=tail(markers$genes, 21), size=2)

knitr::kable(markers)

cell_identity	genes
ISC	Dl
ISC	Smvt
ISC	CG10006
pISC	sna
ISCR2R5	polo
ISC/EB	esg
EB	E(spl)m3-HLH
EB	E(spl)malpha-BFM
EB	E(spl)mbeta-HLH
EB	klu
EE	pros
EE	7B2
EE_subtype	AstA
EE_subtype	AstC
EE_subtype	CCHa1
EE_subtype	CCHa2
EE_subtype	Tk
EE_subtype	NPF
EE_subtype	Orcokinin
EE_subtype	Dh31
aEC	nub
aEC	Myo31DF
mEC	nub
mEC	Myo31DF
pEC	nub
pEC	Myo31DF
aEC	betaTry
pEC	lambdaTry
mEC	lab
mEC	Vha100-4
cardia	Pgant4
copper	PGRP-LA
copper	PGRP-LC
copper	Apoltp
VM	Vn
LFC	Ilp3
LFC	PGRP-SC1a
LFC	PGRP-SC1b
LFC	PGRP-SD
iron	PGRP-SC2
iron	ZIP1(Zip42C.1)

Here we are missing classifying cluster 9, 13, and 16. classying 19 and 20 as ISC cells is pretty weak, IMO, so maybe those too.

clusters = integrated@meta.data$seurat_clusters
integrated_celltype = case_when(
  clusters %in% c(3) ~ "ISC/EB",
  clusters %in% c(12, 14, 15) ~ "EE",
  clusters %in% c(1, 4, 18) ~ "aEC",
  clusters %in% c(6) ~ "iron/copper",
  clusters %in% c(5, 11, 13) ~ "pEC",
  clusters %in% c(8) ~ "mEC",
  clusters %in% c(7) ~ "LFC",
  clusters %in% c(17) ~ "cardia",
  TRUE ~ as.character(clusters))
integrated@meta.data$integrated_celltype = integrated_celltype
Idents(integrated) = integrated@meta.data$integrated_celltype
integratedcelltype.averages = AverageExpression(integrated, return.seurat=TRUE, add.ident="sample")
DoHeatmap(integratedcelltype.averages, features=markers$genes, size=2)

p7 = DimPlot(integrated, group.by="integrated_celltype", label=TRUE)
grid.arrange(p4 + ggtitle("integrated celltypes"), p7 + ggtitle("integrated celltypes"))

knitr::kable(integrated@meta.data %>%
             dplyr::select(seurat_clusters, integrated_celltype) %>%
             unique(), row.names=FALSE)

seurat_clusters	integrated_celltype
7	LFC
5	pEC
2	2
10	10
0	0
11	pEC
6	iron/copper
19	19
8	mEC
21	21
1	aEC
4	aEC
13	pEC
12	EE
3	ISC/EB
14	EE
9	9
16	16
18	aEC
20	20
15	EE
17	cardia

cluster markers (ROC)

Here we will find markers for each of the clusters where we haven’t been able to assign a cell type so we can hopefully figure out what type of cells they are.

Idents(integrated) = integrated@meta.data$seurat_cluster

Cluster 0

dir.create(file.path("results", "markers_roc"))
markersroc_0 = FindMarkers(integrated, ident.1=0, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markersroc_0, file.path("results", "markers_roc", "cluster0-markers_roc.csv"))
knitr::kable(head(markersroc_0, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG12374	0.815	0.7426282	0.630	1.000	0.604
Jon99Ciii	0.790	0.3977139	0.580	0.999	0.669
Bace	0.777	0.4316666	0.554	0.996	0.599
betaTry	0.767	0.5517806	0.534	1.000	0.641
Jon99Cii	0.766	0.3550315	0.532	0.975	0.574
alphaTry	0.759	0.5690638	0.518	1.000	0.790
Jon65Aiii	0.745	0.2787393	0.490	0.975	0.627
CG30025	0.743	0.4160439	0.486	0.973	0.601
pre-rRNA:CR45856	0.257	-1.5008458	0.486	0.133	0.607
Pebp1	0.734	0.3555240	0.468	0.876	0.524
18SrRNA:CR45841	0.273	-1.1015931	0.454	0.158	0.613
18SrRNA-Psi:CR45861	0.273	-1.5898494	0.454	0.091	0.528
18SrRNA:CR41548	0.277	-1.0260107	0.446	0.170	0.623
18SrRNA:CR45838	0.278	-1.0210926	0.444	0.171	0.623
CG12057	0.721	0.6674180	0.442	0.624	0.298
18SrRNA-Psi:CR41602	0.279	-1.5761435	0.442	0.089	0.514
pre-rRNA:CR45846	0.284	-0.9791249	0.432	0.243	0.699
mt:CoI	0.703	0.5018172	0.406	0.984	0.893
mt:CoIII	0.697	0.4695789	0.394	0.961	0.818
His3.3B	0.310	-0.4862776	0.380	0.268	0.682
lncRNA:CR42491	0.311	-0.3863641	0.378	0.189	0.582
14-3-3epsilon	0.316	-0.6482983	0.368	0.169	0.559
crc	0.317	-0.5028142	0.366	0.115	0.479
asRNA:CR45874	0.322	-0.5536836	0.356	0.085	0.421
Hr4	0.323	-0.5593047	0.354	0.142	0.504

Cluster 1

dir.create(file.path("results", "markers_roc"))
markers_roc1 = FindMarkers(integrated, ident.1=1, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc1, file.path("results", "markers_roc", "cluster1-markers_roc.csv"))
knitr::kable(head(markers_roc1, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG6295	0.951	2.4180003	0.902	0.970	0.349
CG12374	0.933	1.7278512	0.866	0.991	0.659
betaTry	0.932	1.4927853	0.864	0.992	0.691
CG30025	0.923	1.5042412	0.846	0.992	0.650
alphaTry	0.919	1.2832729	0.838	0.992	0.819
Jon99Ciii	0.893	1.3578745	0.786	0.989	0.715
Jon99Cii	0.877	1.2946068	0.754	0.978	0.629
CG30031	0.869	1.4176099	0.738	0.954	0.505
gammaTry	0.869	1.4176099	0.738	0.954	0.505
deltaTry	0.861	1.3654433	0.722	0.947	0.502
CG5107	0.853	1.4549836	0.706	0.935	0.488
epsilonTry	0.846	1.4368083	0.692	0.934	0.433
Jon65Aiii	0.824	1.2721982	0.648	0.968	0.675
Bace	0.821	0.9766733	0.642	0.979	0.656
Diedel3	0.809	1.2266437	0.618	0.879	0.467
CG17192	0.808	2.2434487	0.616	0.692	0.104
Jon99Ci	0.800	1.3032560	0.600	0.822	0.325
Amy-p	0.797	0.8226741	0.594	0.848	0.395
Jon25Bi	0.796	2.4973360	0.592	0.712	0.171
Jon65Aiv	0.792	1.0718509	0.584	0.966	0.674
tobi	0.780	1.6524216	0.560	0.749	0.284
Mal-A6	0.763	1.6997284	0.526	0.706	0.250
Acbp5	0.242	-2.2991966	0.516	0.396	0.714
rha	0.747	1.2445401	0.494	0.634	0.189
Mal-A1	0.743	0.9628729	0.486	0.827	0.463

Cluster 2

dir.create(file.path("results", "markers_roc"))
markers_roc2 = FindMarkers(integrated, ident.1=2, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc2, file.path("results", "markers_roc", "cluster2-markers_roc.csv"))
knitr::kable(head(markers_roc2, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG13315	0.781	1.3570550	0.562	0.888	0.628
pre-rRNA:CR45847	0.234	-1.6052927	0.532	0.435	0.771
MtnC	0.758	0.9374882	0.516	0.911	0.673
pre-rRNA:CR45845	0.248	-1.3662956	0.504	0.546	0.838
RpLP1	0.741	0.6907439	0.482	0.939	0.844
RpLP2	0.741	0.5753640	0.482	0.988	0.948
CG34330	0.738	0.8335108	0.476	0.868	0.654
RpL41	0.738	0.6366086	0.476	0.990	0.972
RpS26	0.734	0.6451404	0.468	0.942	0.857
Cam	0.730	0.9412521	0.460	0.872	0.708
RpL37A	0.728	0.6184733	0.456	0.965	0.914
RpL35	0.727	0.6034513	0.454	0.965	0.889
Acbp5	0.724	1.3897104	0.448	0.870	0.672
RpL19	0.724	0.6111552	0.448	0.950	0.862
RpS21	0.723	0.7083381	0.446	0.885	0.773
RpS30	0.720	0.6082303	0.440	0.935	0.857
CG43349	0.719	1.3817990	0.438	0.643	0.296
RpL23	0.717	0.5709788	0.434	0.965	0.920
sta	0.715	0.5939545	0.430	0.949	0.879
RpL15	0.712	0.6074820	0.424	0.948	0.853
RpS27	0.711	0.5947775	0.422	0.954	0.885
pre-rRNA:CR45846	0.290	-1.8422695	0.420	0.354	0.625
RpS28b	0.710	0.6617437	0.420	0.903	0.800
RpL13	0.707	0.6758232	0.414	0.871	0.746
RpL27A	0.706	0.5697464	0.412	0.933	0.863

Cluster 3

markers_roc3 = FindMarkers(integrated, ident.1=3, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc3, file.path("results", "markers_roc", "cluster3-markers_roc.csv"))
knitr::kable(head(markers_roc3, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
lncRNA:CR40469	0.971	2.3022701	0.942	0.999	0.907
lncRNA:CR34335	0.969	2.0554202	0.938	0.996	0.804
alphaTry	0.114	-3.5902805	0.772	0.302	0.877
Vha16-1	0.143	-2.2895249	0.714	0.243	0.837
betaTry	0.149	-3.7497180	0.702	0.116	0.765
Jon65Aiv	0.151	-3.2747248	0.698	0.089	0.748
Jon99Ciii	0.161	-3.7708719	0.678	0.230	0.779
MtnC	0.168	-2.7958055	0.664	0.108	0.738
Bace	0.169	-3.3875524	0.662	0.118	0.728
CG30025	0.172	-3.6942286	0.656	0.121	0.724
CG13315	0.185	-2.2055407	0.630	0.068	0.695
Jon65Aiii	0.185	-2.9903301	0.630	0.205	0.740
CG12374	0.186	-3.4517539	0.628	0.180	0.727
Acbp5	0.188	-2.5010183	0.624	0.159	0.730
RpL41	0.808	0.5931470	0.616	0.992	0.972
Vha13	0.195	-2.0011931	0.610	0.155	0.724
MtnA	0.201	-2.7763653	0.598	0.260	0.744
CG34330	0.202	-1.9493609	0.596	0.176	0.711
CG44008	0.202	-2.0279360	0.596	0.038	0.626
Jon99Cii	0.203	-3.1741526	0.594	0.155	0.699
cib	0.793	1.8690538	0.586	0.716	0.303
bun	0.793	1.4059197	0.586	0.809	0.419
aqz	0.787	0.9852902	0.574	0.869	0.488
Pebp1	0.223	-2.7887394	0.554	0.095	0.638
Tet	0.774	1.9552201	0.548	0.607	0.167

Cluster 4

markers_roc4 = FindMarkers(integrated, ident.1=4, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc4, file.path("results", "markers_roc", "cluster4-markers_roc.csv"))
knitr::kable(head(markers_roc4, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
Bace	0.948	1.6662209	0.896	1.000	0.657
betaTry	0.915	1.3489190	0.830	1.000	0.693
alphaTry	0.914	1.2351138	0.828	1.000	0.820
CG30025	0.908	1.3850072	0.816	1.000	0.652
CG7542	0.902	1.9216820	0.804	0.917	0.277
Diedel3	0.890	1.4456951	0.780	0.986	0.462
Mal-A7	0.874	2.5066108	0.748	0.846	0.268
CG5107	0.873	1.3861378	0.746	0.976	0.488
epsilonTry	0.865	1.3547922	0.730	0.967	0.434
deltaTry	0.861	1.3768449	0.722	0.974	0.504
CG30031	0.861	1.3478686	0.722	0.972	0.507
gammaTry	0.861	1.3478686	0.722	0.972	0.507
Jon65Aiv	0.851	1.0712409	0.702	0.992	0.674
Mal-A1	0.850	1.6094906	0.700	0.919	0.459
Pebp1	0.843	1.2475968	0.686	0.974	0.567
CG8834	0.833	1.2595427	0.666	0.873	0.261
Jon65Aiii	0.825	0.8930227	0.650	0.991	0.676
CG4377	0.812	1.3900148	0.624	0.820	0.302
Cyp4e1	0.808	1.5399858	0.616	0.788	0.314
CG12374	0.777	0.8567660	0.554	0.946	0.665
tobi	0.774	1.5496700	0.548	0.744	0.289
Mal-A8	0.766	1.4995523	0.532	0.663	0.198
MtnA	0.239	-2.7459049	0.522	0.410	0.731
CG4363	0.756	1.3258361	0.512	0.683	0.231
Mal-A6	0.745	1.1122476	0.490	0.677	0.256

Cluster 5

markers_roc5 = FindMarkers(integrated, ident.1=5, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc5, file.path("results", "markers_roc", "cluster5-markers_roc.csv"))
knitr::kable(head(markers_roc5, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
LManVI	0.977	3.2826160	0.954	0.979	0.243
LManV	0.942	2.9089739	0.884	0.915	0.111
CG31343	0.941	2.0399043	0.882	0.972	0.294
ninaD	0.906	2.6390118	0.812	0.871	0.126
LManII	0.904	2.1743168	0.808	0.869	0.141
CG15534	0.887	2.2950064	0.774	0.852	0.135
iotaTry	0.868	1.8477332	0.736	0.836	0.170
Cyp6d5	0.865	1.6516791	0.730	0.891	0.332
Mal-A4	0.859	1.9131132	0.718	0.813	0.161
CG11911	0.855	1.2775764	0.710	0.940	0.346
lambdaTry	0.851	1.5033115	0.702	0.862	0.217
Ance	0.836	1.2035410	0.672	0.882	0.282
LManI	0.834	1.7075755	0.668	0.760	0.109
CG7025	0.834	1.7006117	0.668	0.792	0.132
CG14629	0.831	1.4508569	0.662	0.741	0.135
asRNA:CR45281	0.825	1.1363389	0.650	0.746	0.158
CG9673	0.823	1.0930344	0.646	0.848	0.220
LManIII	0.821	2.2069193	0.642	0.697	0.052
CG15533	0.814	1.4710121	0.628	0.688	0.063
fabp	0.806	0.9226544	0.612	0.929	0.461
CG33127	0.800	1.1443035	0.600	0.836	0.275
Gba1a	0.800	1.1219500	0.600	0.660	0.082
CG11151	0.800	0.9555007	0.600	0.859	0.307
CHMP2B	0.800	0.9366572	0.600	0.704	0.218
CG33306	0.798	1.1217776	0.596	0.718	0.123

Cluster 6

markers_roc6 = FindMarkers(integrated, ident.1=6, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc6, file.path("results", "markers_roc", "cluster6-markers_roc.csv"))
knitr::kable(head(markers_roc6, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
MtnC	0.882	1.7872526	0.764	0.973	0.675
Arc1	0.867	2.4032135	0.734	0.911	0.374
CG43774	0.855	3.1529264	0.710	0.807	0.207
CG34330	0.838	1.4639846	0.676	0.945	0.655
CG5399	0.837	1.6595736	0.674	0.857	0.351
mbl	0.834	2.1404858	0.668	0.850	0.388
MtnD	0.812	2.6718753	0.624	0.775	0.307
MtnA	0.811	1.4680532	0.622	0.959	0.693
CG15423	0.807	1.9412040	0.614	0.737	0.166
CG15422	0.776	1.7420240	0.552	0.730	0.245
Stat92E	0.769	1.7067355	0.538	0.631	0.212
Vha16-1	0.754	0.9391610	0.508	0.952	0.783
MtnB	0.750	1.6650008	0.500	0.708	0.303
Vha13	0.739	0.9406012	0.478	0.900	0.669
Adat1	0.738	1.0285194	0.476	0.640	0.249
thetaTry	0.730	1.9437544	0.460	0.623	0.136
CycG	0.726	1.2622215	0.452	0.812	0.503
MtnE	0.723	1.6779141	0.446	0.642	0.249
pre-rRNA:CR45845	0.279	-1.2473915	0.442	0.621	0.826
RpLP2	0.720	0.6403127	0.440	0.986	0.949
pre-rRNA:CR45847	0.281	-1.3870317	0.438	0.571	0.755
CG8177	0.715	2.2074678	0.430	0.574	0.251
Msp300	0.709	1.2365695	0.418	0.719	0.365
Lk6	0.705	1.2882895	0.410	0.769	0.492
RpL41	0.705	0.5436570	0.410	0.998	0.972

Cluster 7

markers_roc7 = FindMarkers(integrated, ident.1=7, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc7, file.path("results", "markers_roc", "cluster7-markers_roc.csv"))
knitr::kable(head(markers_roc7, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG10472	0.984	2.616301	0.968	1.000	0.405
Jon99Fii	0.934	3.236508	0.868	0.940	0.255
Jon25Biii	0.926	3.256755	0.852	0.960	0.411
CG17571	0.922	2.704814	0.844	0.936	0.257
mag	0.922	2.295322	0.844	0.929	0.172
Acbp5	0.917	1.379220	0.834	0.998	0.670
CG11911	0.904	2.168190	0.808	0.933	0.348
Jon65Ai	0.896	2.885554	0.792	0.858	0.142
Jon99Fi	0.880	3.266696	0.760	0.862	0.228
CG7953	0.861	2.720265	0.722	0.800	0.203
CG16749	0.858	1.553809	0.716	0.894	0.276
CG17633	0.855	2.406228	0.710	0.811	0.222
CG3868	0.849	1.400003	0.698	0.902	0.330
CG8997	0.845	3.178884	0.690	0.801	0.263
Jon65Aii	0.836	3.293820	0.672	0.814	0.210
lncRNA:CR40469	0.170	-2.088592	0.660	0.718	0.925
CG31198	0.829	1.287501	0.658	0.885	0.290
CG7916	0.821	3.257496	0.642	0.750	0.252
Jon99Cii	0.821	1.725367	0.642	0.954	0.642
CG18493	0.820	1.598612	0.640	0.794	0.239
CG31323	0.810	1.136244	0.620	0.836	0.322
Jon99Ciii	0.809	1.736809	0.618	0.965	0.725
CG8093	0.802	1.531794	0.604	0.670	0.092
asRNA:CR45281	0.802	1.095471	0.604	0.709	0.161
CG31233	0.799	1.109342	0.598	0.829	0.263

Cluster 8

markers_roc8 = FindMarkers(integrated, ident.1=8, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc8, file.path("results", "markers_roc", "cluster8-markers_roc.csv"))
knitr::kable(head(markers_roc8, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
Vha16-1	0.965	2.334783	0.930	1.000	0.783
CG5767	0.963	4.182364	0.926	0.949	0.170
Vha13	0.942	1.866896	0.884	0.996	0.667
MtnB	0.940	3.424423	0.880	0.947	0.297
CG30479	0.914	3.367492	0.828	0.898	0.135
Adat1	0.909	2.418973	0.818	0.889	0.242
CG15423	0.907	2.526331	0.814	0.895	0.165
Vha55	0.899	2.306491	0.798	0.958	0.416
CG30480	0.893	2.297150	0.786	0.802	0.056
Vha68-2	0.892	2.124161	0.784	0.931	0.341
CG5399	0.887	2.160550	0.774	0.913	0.354
Vha44	0.886	2.244559	0.772	0.964	0.474
Vha14-1	0.883	1.918986	0.766	0.900	0.299
VhaAC39-1	0.881	2.244777	0.762	0.913	0.328
CG15422	0.868	1.932872	0.736	0.884	0.243
Tsp42Ec	0.865	2.502074	0.730	0.833	0.179
VhaSFD	0.864	2.068645	0.728	0.869	0.244
Vha26	0.849	1.874010	0.698	0.878	0.308
Vha36-1	0.845	1.801660	0.690	0.862	0.304
Argk	0.844	1.768606	0.688	0.909	0.389
VhaAC45	0.838	1.735004	0.676	0.878	0.320
Mpcp2	0.833	1.875708	0.666	0.951	0.558
VhaM9.7-b	0.832	1.618574	0.664	0.842	0.307
CG7430	0.827	2.010335	0.654	0.793	0.292
CAH1	0.824	2.278490	0.648	0.777	0.132

Cluster 9

markers_roc9 = FindMarkers(integrated, ident.1=9, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc9, file.path("results", "markers_roc", "cluster9-markers_roc.csv"))
knitr::kable(head(markers_roc9, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
Amy-p	0.927	3.7546156	0.854	0.926	0.413
Amy-d	0.823	3.4545206	0.646	0.718	0.233
CG11400	0.820	3.3465943	0.640	0.688	0.146
lncRNA:CR40469	0.811	0.6702035	0.622	0.995	0.911
CG4928	0.804	3.3092892	0.608	0.678	0.197
CG3819	0.803	3.3652770	0.606	0.657	0.126
CG15043	0.793	3.0710478	0.586	0.642	0.159
lncRNA:CR34335	0.781	0.7487954	0.562	0.980	0.813
CG6839	0.766	3.6581995	0.532	0.622	0.156
trol	0.755	0.7098187	0.510	0.995	0.990
fmt	0.749	0.6997454	0.498	0.997	0.995
CG14125	0.742	4.9471279	0.484	0.563	0.197
CG13323	0.740	2.4128346	0.480	0.612	0.282
Mal-A1	0.283	-1.8074388	0.434	0.124	0.508
Vha13	0.712	0.6736012	0.424	0.883	0.673
28SrRNA-Psi:CR40741	0.292	-1.5479371	0.416	0.124	0.464
Adhr	0.708	1.6133663	0.416	0.541	0.202
alphaTry	0.295	-1.7875947	0.410	0.761	0.836
CG30025	0.295	-1.8352758	0.410	0.393	0.689
CG15818	0.704	3.2587424	0.408	0.492	0.108
sesB	0.703	0.6889130	0.406	0.893	0.772
Oda	0.699	0.7258382	0.398	0.919	0.830
pAbp	0.695	0.5943689	0.390	0.898	0.761
Vha16-1	0.692	0.2603364	0.384	0.942	0.786
Adh	0.690	1.5147031	0.380	0.503	0.189

Cluster 10

markers_roc10 = FindMarkers(integrated, ident.1=10, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc10, file.path("results", "markers_roc", "cluster10-markers_roc.csv"))
knitr::kable(head(markers_roc10, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
pre-rRNA:CR45846	0.883	1.8543423	0.766	0.967	0.592
pre-rRNA:CR45847	0.877	1.6122745	0.754	0.994	0.737
18SrRNA:CR41548	0.875	1.8092135	0.750	0.917	0.516
18SrRNA:CR45838	0.875	1.8088839	0.750	0.917	0.516
18SrRNA:CR45841	0.869	1.8323186	0.738	0.902	0.506
pre-rRNA:CR45845	0.833	1.5562434	0.666	0.979	0.810
pre-rRNA:CR45856	0.828	1.9342579	0.656	0.842	0.497
18SrRNA-Psi:CR45861	0.807	1.9929792	0.614	0.774	0.425
18SrRNA-Psi:CR41602	0.768	1.9522780	0.536	0.714	0.415
28SrRNA:CR45844	0.715	1.1993126	0.430	0.744	0.607
28SrRNA-Psi:CR45851	0.702	1.3256885	0.404	0.688	0.531
Jon99Cii	0.702	0.5646288	0.404	0.893	0.650
Jon99Ciii	0.681	0.2742590	0.362	0.938	0.731
CG12374	0.678	0.3401385	0.356	0.943	0.678
Bace	0.670	0.2581049	0.340	0.923	0.674
CG30025	0.662	0.2913271	0.324	0.899	0.671
14-3-3epsilon	0.347	-0.6633088	0.306	0.158	0.488
28SrRNA-Psi:CR40741	0.649	1.3577011	0.298	0.565	0.448
28SrRNA-Psi:CR41609	0.642	1.3964064	0.284	0.560	0.459
clos	0.364	-0.3143394	0.272	0.107	0.384
CG33995	0.366	-0.3927518	0.268	0.071	0.341
crc	0.366	-0.3928422	0.268	0.131	0.412
His3.3B	0.366	-0.4291378	0.268	0.304	0.605
CG10472	0.633	0.3546737	0.266	0.601	0.430
CG42322	0.369	-0.2620372	0.262	0.182	0.471

Cluster 11

markers_roc11 = FindMarkers(integrated, ident.1=11, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc11, file.path("results", "markers_roc", "cluster11-markers_roc.csv"))
knitr::kable(head(markers_roc11, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
Gs2	0.960	2.5571275	0.920	0.978	0.369
Acbp3	0.957	2.7989831	0.914	0.978	0.298
CG31343	0.927	1.8074765	0.854	0.984	0.310
zetaTry	0.902	2.2494411	0.804	0.873	0.180
CG13492	0.891	1.8982504	0.782	0.876	0.217
CG9673	0.880	1.7049865	0.760	0.895	0.234
fabp	0.864	1.2838541	0.728	0.962	0.471
CG15254	0.842	2.5365965	0.684	0.796	0.138
MtnA	0.842	1.1704878	0.684	0.987	0.699
CG10116	0.837	1.7645785	0.674	0.812	0.219
Gdh	0.837	1.7188263	0.674	0.768	0.192
CG8774	0.834	2.0324823	0.668	0.768	0.148
Agpat4	0.819	1.5228590	0.638	0.761	0.194
CG32473	0.817	2.6193899	0.634	0.739	0.143
CG5958	0.801	1.7785181	0.602	0.726	0.177
Pepck2	0.795	1.7380777	0.590	0.653	0.096
NAAT1	0.790	1.9975076	0.580	0.669	0.058
dmGlut	0.790	1.4212022	0.580	0.653	0.092
CG10911	0.789	0.9464320	0.578	0.898	0.402
CG32633	0.788	1.2566876	0.576	0.758	0.236
Jon99Ciii	0.215	-3.9366132	0.570	0.455	0.746
Jon65Aiv	0.218	-3.0364897	0.564	0.366	0.708
CG4653	0.781	1.4001194	0.562	0.736	0.220
Acbp5	0.779	0.7779398	0.558	0.978	0.679
Jon65Aiii	0.221	-3.0295947	0.558	0.360	0.710

Cluster 12

markers_roc12 = FindMarkers(integrated, ident.1=12, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc12, file.path("results", "markers_roc", "cluster12-markers_roc.csv"))
knitr::kable(head(markers_roc12, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
IA-2	0.964	2.647187	0.928	0.993	0.403
nrv3	0.916	2.247964	0.832	0.909	0.242
7B2	0.913	2.096422	0.826	0.915	0.231
Phm	0.882	1.793405	0.764	0.866	0.234
Pal2	0.852	2.079188	0.704	0.782	0.147
pros	0.844	2.015971	0.688	0.769	0.181
w	0.841	2.291362	0.682	0.759	0.206
CG30183	0.838	1.964119	0.676	0.736	0.130
alphaTry	0.172	-3.018572	0.656	0.440	0.845
Jon99Ciii	0.179	-3.040954	0.642	0.212	0.753
Hsp83	0.813	1.332282	0.626	0.850	0.374
Jon65Aiv	0.195	-2.829406	0.610	0.186	0.714
betaTry	0.195	-3.142426	0.610	0.215	0.731
heph	0.804	1.849458	0.608	0.664	0.104
Galphao	0.803	1.910813	0.606	0.691	0.175
Jon65Aiii	0.200	-2.671267	0.600	0.199	0.715
RpL8	0.201	-1.069680	0.598	0.560	0.903
Gbeta13F	0.798	1.369896	0.596	0.779	0.335
Bace	0.202	-2.942945	0.596	0.160	0.698
RpL35	0.203	-1.058824	0.594	0.564	0.904
RpL37A	0.207	-1.029812	0.586	0.664	0.925
CG12374	0.207	-3.043524	0.586	0.189	0.701
RpL17	0.211	-1.002210	0.578	0.560	0.895
Jon99Cii	0.216	-2.491872	0.568	0.147	0.673
CG30025	0.216	-3.165431	0.568	0.238	0.691

Cluster 13

markers_roc13 = FindMarkers(integrated, ident.1=13, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc13, file.path("results", "markers_roc", "cluster13-markers_roc.csv"))
knitr::kable(head(markers_roc13, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG33926	0.924	2.1060157	0.848	0.950	0.228
CG10912	0.911	2.2182013	0.822	0.907	0.203
CG10911	0.891	1.5862509	0.782	0.963	0.401
CG32368	0.883	2.1692952	0.766	0.854	0.168
Mur29B	0.855	2.5662759	0.710	0.804	0.184
CG42825	0.844	2.2269327	0.688	0.761	0.116
CG14499	0.843	3.0172332	0.686	0.764	0.075
lectin-37Da	0.839	2.7082242	0.678	0.757	0.081
CG31086	0.823	1.5386260	0.646	0.804	0.311
CG31323	0.821	1.7430702	0.642	0.837	0.335
CG18493	0.794	1.4654841	0.588	0.754	0.253
RpL3	0.213	-0.9937845	0.574	0.718	0.896
eEF1alpha1	0.218	-0.8581779	0.564	0.887	0.970
Npc2e	0.778	2.1977119	0.556	0.691	0.154
CG15255	0.778	1.3555371	0.556	0.771	0.267
asRNA:CR44192	0.778	1.3426450	0.556	0.635	0.150
Fer2LCH	0.772	0.9810805	0.544	0.897	0.479
RpS29	0.238	-0.8599220	0.524	0.837	0.941
lectin-37Db	0.757	1.3335771	0.514	0.688	0.226
CG9568	0.754	1.6351462	0.508	0.688	0.224
CG44008	0.754	0.6956464	0.508	0.924	0.571
RpS2	0.247	-0.9726861	0.506	0.731	0.873
sta	0.249	-0.9397450	0.502	0.721	0.889
RpL14	0.250	-0.9298185	0.500	0.638	0.842
RpL38	0.254	-0.9220046	0.492	0.585	0.843

Cluster 14

markers_roc14 = FindMarkers(integrated, ident.1=14, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc14, file.path("results", "markers_roc", "cluster14-markers_roc.csv"))
knitr::kable(head(markers_roc14, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
NPF	0.989	5.147902	0.978	0.983	0.290
IA-2	0.986	3.213505	0.972	1.000	0.403
7B2	0.933	3.044842	0.866	0.924	0.231
Phm	0.913	2.665425	0.826	0.890	0.234
Tk	0.898	2.448163	0.796	0.857	0.210
svr	0.896	2.567417	0.792	0.867	0.243
chrb	0.884	3.500329	0.768	0.841	0.281
nrv3	0.867	2.251026	0.734	0.827	0.245
CG30183	0.852	2.181997	0.704	0.764	0.129
CG46385	0.850	2.210807	0.700	0.831	0.325
esg	0.840	2.222163	0.680	0.761	0.184
unc-13-4A	0.834	1.715761	0.668	0.757	0.169
Hsp23	0.827	2.841930	0.654	0.781	0.318
Ldh	0.827	2.230674	0.654	0.728	0.177
alphaTry	0.177	-2.887524	0.646	0.482	0.844
cpo	0.822	2.446349	0.644	0.761	0.275
Pal2	0.822	2.156723	0.644	0.748	0.149
Hsp22	0.822	1.852289	0.644	0.777	0.341
mbl	0.821	1.258315	0.642	0.867	0.399
Hsp83	0.820	1.744655	0.640	0.847	0.374
Unr	0.815	1.970741	0.630	0.804	0.377
Hsp26	0.809	1.839010	0.618	0.847	0.432
Ih	0.807	2.434480	0.614	0.661	0.097
shep	0.807	1.896214	0.614	0.728	0.241
CG30025	0.195	-2.890765	0.610	0.199	0.692

Cluster 15

markers_roc15 = FindMarkers(integrated, ident.1=15, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc15, file.path("results", "markers_roc", "cluster15-markers_roc.csv"))
knitr::kable(head(markers_roc15, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
Orcokinin	1.000	5.327171	1.000	1.000	0.292
AstC	0.987	3.333166	0.974	0.996	0.293
svr	0.965	3.009325	0.930	0.979	0.244
unc-13-4A	0.957	2.946501	0.914	0.963	0.167
nrv3	0.939	2.642404	0.878	0.954	0.246
IA-2	0.936	2.182955	0.872	0.992	0.407
7B2	0.927	2.126905	0.854	0.959	0.235
Phm	0.923	2.484770	0.846	0.934	0.237
tap	0.917	1.980146	0.834	0.851	0.030
heph	0.898	2.319493	0.796	0.846	0.103
CG30183	0.898	1.741654	0.796	0.880	0.130
pros	0.896	2.421778	0.792	0.871	0.183
CG14989	0.884	3.799407	0.768	0.871	0.190
Ca-alpha1T	0.880	2.284545	0.760	0.788	0.055
slo	0.880	1.877297	0.760	0.780	0.045
Rgk3	0.879	2.170133	0.758	0.768	0.022
jus	0.873	1.609255	0.746	0.768	0.051
Pal2	0.871	2.284061	0.742	0.846	0.150
epsilonTry	0.129	-2.541908	0.742	0.046	0.484
lncRNA:CR31451	0.870	1.543013	0.740	0.780	0.062
fkh	0.867	1.852307	0.734	0.871	0.279
CG44247	0.867	1.801113	0.734	0.780	0.060
shep	0.864	1.872548	0.728	0.876	0.241
CG46385	0.859	1.946008	0.718	0.905	0.327
fru	0.859	1.606407	0.718	0.780	0.144

Cluster 16

markers_roc16 = FindMarkers(integrated, ident.1=16, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc16, file.path("results", "markers_roc", "cluster16-markers_roc.csv"))
knitr::kable(head(markers_roc16, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
EbpIII	0.948	4.834704	0.896	0.920	0.223
Acbp2	0.914	2.549511	0.828	0.903	0.364
lncRNA:CR34335	0.905	1.580156	0.810	1.000	0.815
Adhr	0.883	3.095023	0.766	0.819	0.200
Adh	0.861	3.083127	0.722	0.776	0.187
lncRNA:CR40469	0.858	1.202531	0.716	1.000	0.912
to	0.853	3.587573	0.706	0.726	0.064
spidey	0.841	3.329530	0.682	0.734	0.216
CG30197	0.825	3.532219	0.650	0.684	0.094
CG31523	0.824	2.969051	0.648	0.658	0.148
CG31522	0.807	3.284502	0.614	0.688	0.226
CG1124	0.786	3.623984	0.572	0.620	0.110
Msr-110	0.779	2.666640	0.558	0.612	0.116
CG6770	0.768	1.183859	0.536	0.819	0.557
ple	0.763	2.881654	0.526	0.553	0.030
wat	0.756	3.231272	0.512	0.532	0.046
CG10237	0.750	2.613266	0.500	0.536	0.048
Moe	0.748	1.735258	0.496	0.591	0.240
emp	0.744	2.438345	0.488	0.523	0.089
vir-1	0.738	2.554232	0.476	0.494	0.033
Hacd1	0.737	2.129050	0.474	0.485	0.111
CG8306	0.734	2.509803	0.468	0.494	0.060
ATPCL	0.730	2.223677	0.460	0.527	0.170
Inos	0.729	2.157946	0.458	0.502	0.099
Taldo	0.727	1.879034	0.454	0.523	0.197

Cluster 17

markers_roc17 = FindMarkers(integrated, ident.1=17, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc17, file.path("results", "markers_roc", "cluster17-markers_roc.csv"))
knitr::kable(head(markers_roc17, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG14645	0.999	5.166353	0.998	1.000	0.282
Skp2	0.998	5.300042	0.996	1.000	0.389
CG3906	0.992	5.087155	0.984	0.991	0.217
CG34324	0.990	5.281636	0.980	0.987	0.277
CG34220	0.988	5.093672	0.976	0.987	0.420
CG4783	0.971	3.204050	0.942	0.956	0.160
Muc68D	0.958	5.030528	0.916	0.933	0.207
CG11672	0.932	3.842805	0.864	0.871	0.037
Idgf4	0.874	2.455162	0.748	0.782	0.072
Pgant4	0.828	2.450677	0.656	0.662	0.014
alphaTry	0.186	-2.869711	0.628	0.440	0.842
CG30025	0.186	-3.163239	0.628	0.116	0.690
Jon99Ciii	0.201	-2.909744	0.598	0.227	0.749
CG34330	0.214	-2.012285	0.572	0.133	0.682
Jon65Aiv	0.216	-2.543440	0.568	0.204	0.709
Jon65Aiii	0.217	-2.454261	0.566	0.209	0.710
betaTry	0.217	-2.727593	0.566	0.236	0.726
Jon99Cii	0.229	-2.470175	0.542	0.151	0.669
Bace	0.232	-2.477221	0.536	0.209	0.693
Vha16-1	0.239	-1.611035	0.522	0.364	0.801
MtnC	0.239	-2.133976	0.522	0.280	0.700
sesB	0.240	-1.164468	0.520	0.320	0.787
Act5C	0.243	-1.032861	0.514	0.604	0.900
CG12374	0.245	-2.491602	0.510	0.249	0.696
CG30026	0.754	1.759105	0.508	0.511	0.003

Cluster 18

markers_roc18 = FindMarkers(integrated, ident.1=18, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc18, file.path("results", "markers_roc", "cluster18-markers_roc.csv"))
knitr::kable(head(markers_roc18, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
Npc2f	0.986	3.6101816	0.972	0.988	0.150
Jon65Aiv	0.986	2.1546865	0.972	1.000	0.694
thetaTry	0.980	3.0983381	0.960	0.988	0.149
yip7	0.979	2.3389780	0.958	1.000	0.539
Jon65Aiii	0.979	2.0314714	0.958	1.000	0.695
Peritrophin-15a	0.938	2.6028515	0.876	0.950	0.184
Bace	0.926	1.8774625	0.852	0.994	0.678
CG30031	0.915	1.3828164	0.830	1.000	0.535
gammaTry	0.915	1.3828164	0.830	1.000	0.535
deltaTry	0.890	1.3582514	0.780	1.000	0.532
CG4734	0.885	1.8411515	0.770	0.826	0.132
CG18404	0.879	4.0870292	0.758	0.845	0.161
alphaTry	0.872	1.0091274	0.744	1.000	0.831
Pebp1	0.867	1.5259964	0.734	0.944	0.592
CG9682	0.851	2.5655588	0.702	0.727	0.027
CG30025	0.842	1.0353096	0.684	1.000	0.673
LysB	0.835	1.3202083	0.670	0.826	0.197
CG4830	0.827	2.6994819	0.654	0.789	0.017
lncRNA:CR40469	0.181	-2.0822384	0.638	0.826	0.915
CG4563	0.816	2.4116003	0.632	0.770	0.061
CG5107	0.799	0.8569231	0.598	0.963	0.518
betaTry	0.788	0.8541648	0.576	1.000	0.711
CG34040	0.783	0.5422801	0.566	0.640	0.082
CG4377	0.779	0.9930001	0.558	0.857	0.333
epsilonTry	0.778	0.7924996	0.556	0.957	0.467

Cluster 19

markers_roc19 = FindMarkers(integrated, ident.1=19, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc19, file.path("results", "markers_roc", "cluster19-markers_roc.csv"))
knitr::kable(head(markers_roc19, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG6409	0.980	4.8036985	0.960	0.966	0.145
CG42656	0.951	4.4368015	0.902	0.910	0.082
IM33	0.908	3.9650892	0.816	0.831	0.081
Vha13	0.906	1.5702548	0.812	0.978	0.679
lncRNA:CR40469	0.902	1.3818689	0.804	1.000	0.913
Vha16-1	0.892	1.1695161	0.784	0.978	0.791
sesB	0.880	1.2228950	0.760	0.978	0.775
Vha55	0.845	1.5750791	0.690	0.876	0.436
lncRNA:CR34335	0.842	1.0495865	0.684	0.989	0.817
ATPsynC	0.817	1.1686827	0.634	0.910	0.659
Vha44	0.815	1.3559872	0.630	0.843	0.492
Gdh	0.811	2.5818285	0.622	0.685	0.205
Vha26	0.808	1.6921151	0.616	0.764	0.328
Vha68-2	0.805	1.5239685	0.610	0.775	0.362
porin	0.799	1.5868938	0.598	0.809	0.503
tws	0.797	2.6627324	0.594	0.640	0.170
blw	0.790	1.4953967	0.580	0.820	0.556
Ggamma30A	0.778	2.4257360	0.556	0.584	0.061
CG1143	0.776	2.8464808	0.552	0.584	0.052
Irk2	0.775	2.8819032	0.550	0.584	0.057
Msr-110	0.764	2.2743860	0.528	0.596	0.123
pAbp	0.760	0.7985467	0.520	0.910	0.765
Argk	0.753	1.2657596	0.506	0.730	0.408
Mpcp2	0.753	1.1475929	0.506	0.798	0.572
Dic1	0.752	1.8342125	0.504	0.562	0.178

Cluster 20

markers_roc20 = FindMarkers(integrated, ident.1=20, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc20, file.path("results", "markers_roc", "cluster20-markers_roc.csv"))
knitr::kable(head(markers_roc20, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
CG13285	0.980	3.8504561	0.960	0.963	0.020
Msr-110	0.970	3.1788920	0.940	0.963	0.123
CG44013	0.887	2.9005776	0.774	0.778	0.020
lncRNA:CR34335	0.886	1.2155838	0.772	1.000	0.818
Idgf4	0.885	2.3408393	0.770	0.815	0.083
Idgf6	0.881	2.6238014	0.762	0.778	0.043
lncRNA:CR40469	0.876	0.9434807	0.752	1.000	0.913
CG5065	0.852	2.2587751	0.704	0.741	0.095
CG43394	0.837	2.6414936	0.674	0.685	0.028
CG10096	0.829	2.8104611	0.658	0.667	0.012
mt:lrRNA	0.819	0.3111689	0.638	1.000	1.000
Vha16-1	0.186	-2.1373285	0.628	0.241	0.795
MtnC	0.201	-2.6273170	0.598	0.111	0.694
Act5C	0.202	-1.2459382	0.596	0.574	0.896
Orcokinin	0.794	0.5212523	0.588	0.741	0.306
alphaTry	0.208	-2.5242600	0.584	0.444	0.836
Jon65Aiii	0.212	-2.6360940	0.576	0.148	0.702
Jon99Cii	0.227	-2.9229030	0.546	0.093	0.661
GstD11	0.768	2.0589475	0.536	0.537	0.001
Jon65Aiv	0.232	-2.4779429	0.536	0.204	0.701
Acbp5	0.234	-1.7272450	0.532	0.148	0.690
Cys	0.765	1.8452262	0.530	0.611	0.177
CG30025	0.240	-2.5843813	0.520	0.204	0.681
Vha13	0.242	-1.7401562	0.516	0.167	0.684
Atpalpha	0.754	1.4416261	0.508	0.667	0.310

Cluster 21

markers_roc21 = FindMarkers(integrated, ident.1=21, test.use="roc") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(-power)
readr::write_csv(markers_roc21, file.path("results", "markers_roc", "cluster21-markers_roc.csv"))
knitr::kable(head(markers_roc21, 25))

gene	myAUC	avg_diff	power	pct.1	pct.2
MtnA	0.933	1.8253924	0.866	1.000	0.706
Alp4	0.810	4.0091710	0.620	0.639	0.055
Vha13	0.801	1.0715504	0.602	0.889	0.681
CG6726	0.784	3.1950727	0.568	0.611	0.109
Mur18B	0.773	3.7634567	0.546	0.583	0.080
Alp2	0.771	3.2290870	0.542	0.556	0.029
sesB	0.756	0.7051319	0.512	0.944	0.776
CG3168	0.746	3.7752656	0.492	0.528	0.073
Vha26	0.745	1.0301148	0.490	0.694	0.330
alphaTry	0.264	-1.4814880	0.472	0.556	0.835
HDAC6	0.729	2.5162694	0.458	0.500	0.161
Mal-A1	0.272	-2.6635592	0.456	0.028	0.495
CG18324	0.723	1.7456577	0.446	0.444	0.126
UGP	0.721	2.2285553	0.442	0.500	0.160
Vha16-1	0.721	0.6729296	0.442	0.889	0.792
CG30025	0.286	-1.8776482	0.428	0.333	0.679
smp-30	0.707	2.1958556	0.414	0.472	0.166
Vha44	0.703	0.8803877	0.406	0.694	0.494
CG31288	0.700	2.0329887	0.400	0.444	0.172
CG30031	0.301	-1.3683171	0.398	0.167	0.543
gammaTry	0.301	-1.3683171	0.398	0.167	0.543
CG13868	0.698	1.8833634	0.396	0.417	0.137
deltaTry	0.303	-1.3374017	0.394	0.167	0.540
CG5107	0.303	-2.0925566	0.394	0.139	0.526
Vha55	0.695	0.8688583	0.390	0.667	0.439

Specific markers-roc plot for cluster 21.

DimPlot(integrated, label=TRUE)

FeaturePlot(integrated, features=c("MtnA", "Vha13", "Smvt", "Oatp58Dc"))

cells21 = rownames(subset(integrated@meta.data, seurat_clusters == 21))
FeaturePlot(integrated[,cells21], features=c("MtnA", "Vha13", "Smvt", "Oatp58Dc")) +
  ggtitle("cluster 21 only")

cluster markers (MAST)

Here we will find markers for each of the clusters where we haven’t been able to assign a cell type so we can hopefully figure out what type of cells they are.

Cluster 0

dir.create(file.path("results", "markers_MAST"))
markers_MAST0 = FindMarkers(integrated, ident.1=0, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST0, file.path("results", "markers_MAST", "cluster0-markers_MAST.csv"))
knitr::kable(head(markers_MAST0, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CG12374	0	0.7426282	1.000	0.604	0
alphaTry	0	0.5690638	1.000	0.790	0
betaTry	0	0.5517806	1.000	0.641	0
Bace	0	0.4316666	0.996	0.599	0
CG30025	0	0.4160439	0.973	0.601	0
Jon99Ciii	0	0.3977139	0.999	0.669	0
Jon99Cii	0	0.3550315	0.975	0.574	0
Jon65Aiii	0	0.2787393	0.975	0.627	0
pre-rRNA:CR45846	0	-0.9791249	0.243	0.699	0
18SrRNA:CR45838	0	-1.0210926	0.171	0.623	0
18SrRNA:CR41548	0	-1.0260107	0.170	0.623	0
18SrRNA:CR45841	0	-1.1015931	0.158	0.613	0
pre-rRNA:CR45856	0	-1.5008458	0.133	0.607	0
18SrRNA-Psi:CR41602	0	-1.5761435	0.089	0.514	0
18SrRNA-Psi:CR45861	0	-1.5898494	0.091	0.528	0
CG12057	0	0.6674180	0.624	0.298	0
His3.3B	0	-0.4862776	0.268	0.682	0
Pebp1	0	0.3555240	0.876	0.524	0
ps	0	-0.3731687	0.233	0.625	0
lncRNA:CR42491	0	-0.3863641	0.189	0.582	0
14-3-3epsilon	0	-0.6482983	0.169	0.559	0
cib	0	-0.5655328	0.066	0.405	0
CycG	0	-0.8424565	0.207	0.601	0
porin	0	-0.3450815	0.217	0.582	0
28SrRNA-Psi:CR45862	0	-0.2617602	0.030	0.339	0

Cluster 1

dir.create(file.path("results", "markers_MAST"))
markers_MAST1 = FindMarkers(integrated, ident.1=1, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST1, file.path("results", "markers_MAST", "cluster1-markers_MAST.csv"))
knitr::kable(head(markers_MAST1, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Jon25Bi	0	2.4973360	0.712	0.171	0
CG6295	0	2.4180003	0.970	0.349	0
CG17192	0	2.2434487	0.692	0.104	0
CG34026	0	1.8928975	0.535	0.076	0
CG12374	0	1.7278512	0.991	0.659	0
CG30025	0	1.5042412	0.992	0.650	0
betaTry	0	1.4927853	0.992	0.691	0
CG30031	0	1.4176099	0.954	0.505	0
gammaTry	0	1.4176099	0.954	0.505	0
deltaTry	0	1.3654433	0.947	0.502	0
CG5107	0	1.4549836	0.935	0.488	0
Jon99Cii	0	1.2946068	0.978	0.629	0
epsilonTry	0	1.4368083	0.934	0.433	0
Jon99Ciii	0	1.3578745	0.989	0.715	0
Jon25Bii	0	2.2726902	0.610	0.146	0
alphaTry	0	1.2832729	0.992	0.819	0
Jon99Ci	0	1.3032560	0.822	0.325	0
Mal-A6	0	1.6997284	0.706	0.250	0
Jon65Aiii	0	1.2721982	0.968	0.675	0
Bace	0	0.9766733	0.979	0.656	0
Diedel3	0	1.2266437	0.879	0.467	0
MtnA	0	-2.4875220	0.519	0.724	0
Amy-p	0	0.8226741	0.848	0.395	0
tobi	0	1.6524216	0.749	0.284	0
Acbp5	0	-2.2991966	0.396	0.714	0

Cluster 2

dir.create(file.path("results", "markers_MAST"))
markers_MAST2 = FindMarkers(integrated, ident.1=2, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST2, file.path("results", "markers_MAST", "cluster2-markers_MAST.csv"))
knitr::kable(head(markers_MAST2, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
pre-rRNA:CR45847	0	-1.6052927	0.435	0.771	0
CG13315	0	1.3570550	0.888	0.628	0
pre-rRNA:CR45845	0	-1.3662956	0.546	0.838	0
CG43349	0	1.3817990	0.643	0.296	0
pre-rRNA:CR45846	0	-1.8422695	0.354	0.625	0
MtnC	0	0.9374882	0.911	0.673	0
RpLP1	0	0.6907439	0.939	0.844	0
RpS26	0	0.6451404	0.942	0.857	0
Acbp5	0	1.3897104	0.870	0.672	0
RpS21	0	0.7083381	0.885	0.773	0
Cam	0	0.9412521	0.872	0.708	0
RpLP2	0	0.5753640	0.988	0.948	0
CG34330	0	0.8335108	0.868	0.654	0
28SrRNA:CR45844	0	-1.1858404	0.391	0.629	0
RpL41	0	0.6366086	0.990	0.972	0
RpL37A	0	0.6184733	0.965	0.914	0
RpL19	0	0.6111552	0.950	0.862	0
RpS30	0	0.6082303	0.935	0.857	0
RpS27	0	0.5947775	0.954	0.885	0
RpL13	0	0.6758232	0.871	0.746	0
RpL35	0	0.6034513	0.965	0.889	0
RpL23	0	0.5709788	0.965	0.920	0
RpS28b	0	0.6617437	0.903	0.800	0
18SrRNA:CR45838	0	-1.4476648	0.340	0.544	0
18SrRNA:CR41548	0	-1.4466951	0.340	0.544	0

Cluster 3

markers_MAST3 = FindMarkers(integrated, ident.1=3, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST3, file.path("results", "markers_MAST", "cluster3-markers_MAST.csv"))
knitr::kable(head(markers_MAST3, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
lncRNA:CR40469	0	2.302270	0.999	0.907	0
lncRNA:CR34335	0	2.055420	0.996	0.804	0
Tet	0	1.955220	0.607	0.167	0
sty	0	1.890713	0.616	0.261	0
cib	0	1.869054	0.716	0.303	0
E(spl)mbeta-HLH	0	1.832528	0.498	0.121	0
Df31	0	1.714096	0.662	0.234	0
E(spl)malpha-BFM	0	1.699416	0.438	0.115	0
hdc	0	1.682541	0.539	0.210	0
alphaTry	0	-3.590280	0.302	0.877	0
Jon65Aiv	0	-3.274725	0.089	0.748	0
esg	0	1.623873	0.627	0.165	0
betaTry	0	-3.749718	0.116	0.765	0
lncRNA:CR44898	0	1.181112	0.417	0.093	0
MtnC	0	-2.795806	0.108	0.738	0
mew	0	1.489866	0.524	0.196	0
robo2	0	1.338593	0.374	0.084	0
Vha16-1	0	-2.289525	0.243	0.837	0
CG30025	0	-3.694229	0.121	0.724	0
Jon99Ciii	0	-3.770872	0.230	0.779	0
Bace	0	-3.387552	0.118	0.728	0
CG13315	0	-2.205541	0.068	0.695	0
N	0	1.341854	0.367	0.102	0
Jon99Cii	0	-3.174153	0.155	0.699	0
Jon65Aiii	0	-2.990330	0.205	0.740	0

Cluster 4

markers_MAST4 = FindMarkers(integrated, ident.1=4, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST4, file.path("results", "markers_MAST", "cluster4-markers_MAST.csv"))
knitr::kable(head(markers_MAST4, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Mal-A7	0	2.5066108	0.846	0.268	0
CG7542	0	1.9216820	0.917	0.277	0
Bace	0	1.6662209	1.000	0.657	0
Diedel3	0	1.4456951	0.986	0.462	0
CG5107	0	1.3861378	0.976	0.488	0
epsilonTry	0	1.3547922	0.967	0.434	0
deltaTry	0	1.3768449	0.974	0.504	0
CG30031	0	1.3478686	0.972	0.507	0
gammaTry	0	1.3478686	0.972	0.507	0
CG30025	0	1.3850072	1.000	0.652	0
betaTry	0	1.3489190	1.000	0.693	0
CG8834	0	1.2595427	0.873	0.261	0
Mal-A1	0	1.6094906	0.919	0.459	0
Cyp4e1	0	1.5399858	0.788	0.314	0
alphaTry	0	1.2351138	1.000	0.820	0
Pebp1	0	1.2475968	0.974	0.567	0
CG4377	0	1.3900148	0.820	0.302	0
Jon65Aiv	0	1.0712409	0.992	0.674	0
Mal-A8	0	1.4995523	0.663	0.198	0
CG4363	0	1.3258361	0.683	0.231	0
Jon65Aiii	0	0.8930227	0.991	0.676	0
tobi	0	1.5496700	0.744	0.289	0
MtnA	0	-2.7459049	0.410	0.731	0
LysD	0	1.1315302	0.605	0.207	0
CG12374	0	0.8567660	0.946	0.665	0

Cluster 5

markers_MAST5 = FindMarkers(integrated, ident.1=5, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST5, file.path("results", "markers_MAST", "cluster5-markers_MAST.csv"))
knitr::kable(head(markers_MAST5, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
LManVI	0	3.2826160	0.979	0.243	0
LManV	0	2.9089739	0.915	0.111	0
ninaD	0	2.6390118	0.871	0.126	0
CG15534	0	2.2950064	0.852	0.135	0
LManIII	0	2.2069193	0.697	0.052	0
LManII	0	2.1743168	0.869	0.141	0
CG31343	0	2.0399043	0.972	0.294	0
Mal-A4	0	1.9131132	0.813	0.161	0
iotaTry	0	1.8477332	0.836	0.170	0
LManI	0	1.7075755	0.760	0.109	0
CG7025	0	1.7006117	0.792	0.132	0
Jon99Cii	0	-2.6760089	0.608	0.660	0
CG14629	0	1.4508569	0.741	0.135	0
yip7	0	-2.6590889	0.577	0.544	0
Jon65Aiv	0	-2.9567791	0.605	0.704	0
CG15533	0	1.4710121	0.688	0.063	0
Pebp1	0	-2.5958549	0.586	0.598	0
Cyp6d5	0	1.6516791	0.891	0.332	0
Jon65Aiii	0	-2.9732678	0.612	0.705	0
lambdaTry	0	1.5033115	0.862	0.217	0
Npc2d	0	1.7057162	0.684	0.113	0
Rrp46	0	0.2793078	0.541	0.033	0
CG11911	0	1.2775764	0.940	0.346	0
asRNA:CR45281	0	1.1363389	0.746	0.158	0
Jon99Ciii	0	-3.2440617	0.649	0.743	0

Cluster 6

markers_MAST6 = FindMarkers(integrated, ident.1=6, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST6, file.path("results", "markers_MAST", "cluster6-markers_MAST.csv"))
knitr::kable(head(markers_MAST6, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CG43774	0	3.1529264	0.807	0.207	0
MtnD	0	2.6718753	0.775	0.307	0
Arc1	0	2.4032135	0.911	0.374	0
MtnC	0	1.7872526	0.973	0.675	0
CG15423	0	1.9412040	0.737	0.166	0
mbl	0	2.1404858	0.850	0.388	0
thetaTry	0	1.9437544	0.623	0.136	0
CG34330	0	1.4639846	0.945	0.655	0
CG5773	0	1.8438525	0.352	0.028	0
CG5399	0	1.6595736	0.857	0.351	0
CG34301	0	2.0971596	0.504	0.156	0
Gagr	0	0.5444069	0.415	0.024	0
CG15422	0	1.7420240	0.730	0.245	0
CG8177	0	2.2074678	0.574	0.251	0
CG15127	0	1.7561668	0.354	0.076	0
Stat92E	0	1.7067355	0.631	0.212	0
MtnB	0	1.6650008	0.708	0.303	0
CG5770	0	1.6160867	0.315	0.023	0
MtnA	0	1.4680532	0.959	0.693	0
pre-rRNA:CR45847	0	-1.3870317	0.571	0.755	0
MtnE	0	1.6779141	0.642	0.249	0
28SrRNA:CR45844	0	-1.2991613	0.504	0.617	0
Tsp42Ec	0	1.5186541	0.562	0.187	0
CG9672	0	1.7261596	0.306	0.112	0
Adat1	0	1.0285194	0.640	0.249	0

Cluster 7

markers_MAST7 = FindMarkers(integrated, ident.1=7, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST7, file.path("results", "markers_MAST", "cluster7-markers_MAST.csv"))
knitr::kable(head(markers_MAST7, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Jon65Aii	0	3.293820	0.814	0.210	0
Jon99Fi	0	3.266696	0.862	0.228	0
CG7916	0	3.257496	0.750	0.252	0
Jon25Biii	0	3.256755	0.960	0.411	0
Jon99Fii	0	3.236508	0.940	0.255	0
Jon65Ai	0	2.885554	0.858	0.142	0
CG7953	0	2.720265	0.800	0.203	0
CG17571	0	2.704814	0.936	0.257	0
CG10472	0	2.616301	1.000	0.405	0
mag	0	2.295322	0.929	0.172	0
CG3106	0	1.981872	0.638	0.041	0
CG8997	0	3.178884	0.801	0.263	0
Ag5r2	0	1.913701	0.658	0.116	0
CG11911	0	2.168190	0.933	0.348	0
CG17633	0	2.406228	0.811	0.222	0
PGRP-SC1b	0	0.844151	0.454	0.007	0
CG16749	0	1.553809	0.894	0.276	0
Acbp5	0	1.379220	0.998	0.670	0
CG8093	0	1.531794	0.670	0.092	0
Jon66Ci	0	2.028641	0.552	0.041	0
Jon44E	0	2.093519	0.517	0.065	0
CG3868	0	1.400003	0.902	0.330	0
Bace	0	-3.270406	0.466	0.694	0
CG31198	0	1.287501	0.885	0.290	0
asRNA:CR45281	0	1.095471	0.709	0.161	0

Cluster 8

markers_MAST8 = FindMarkers(integrated, ident.1=8, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST8, file.path("results", "markers_MAST", "cluster8-markers_MAST.csv"))
knitr::kable(head(markers_MAST8, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CG5767	0	4.1823644	0.949	0.170	0
MtnB	0	3.4244227	0.947	0.297	0
CG30479	0	3.3674916	0.898	0.135	0
CG6277	0	3.3655225	0.646	0.112	0
CG8661	0	3.3205974	0.757	0.145	0
Vha100-4	0	2.9384083	0.664	0.047	0
CG31446	0	2.6010173	0.653	0.064	0
CG15423	0	2.5263314	0.895	0.165	0
Tsp42Ec	0	2.5020736	0.833	0.179	0
Adat1	0	2.4189735	0.889	0.242	0
Vha16-1	0	2.3347834	1.000	0.783	0
CG30480	0	2.2971498	0.802	0.056	0
Vha13	0	1.8668964	0.996	0.667	0
CAH1	0	2.2784895	0.777	0.132	0
CG17109	0	2.0909091	0.644	0.130	0
Vha55	0	2.3064912	0.958	0.416	0
CG5399	0	2.1605505	0.913	0.354	0
CG11192	0	0.3031397	0.474	0.005	0
CG31087	0	2.0721848	0.751	0.208	0
Vha14-1	0	1.9189858	0.900	0.299	0
Vha44	0	2.2445595	0.964	0.474	0
Vha68-2	0	2.1241612	0.931	0.341	0
CG15422	0	1.9328725	0.884	0.243	0
VhaAC39-1	0	2.2447769	0.913	0.328	0
CG6271	0	1.1021336	0.428	0.005	0

Cluster 9

markers_MAST9 = FindMarkers(integrated, ident.1=9, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST9, file.path("results", "markers_MAST", "cluster9-markers_MAST.csv"))
knitr::kable(head(markers_MAST9, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Amy-p	0	3.7546156	0.926	0.413	0
CG6839	0	3.6581995	0.622	0.156	0
CG3819	0	3.3652770	0.657	0.126	0
CG11400	0	3.3465943	0.688	0.146	0
CG4928	0	3.3092892	0.678	0.197	0
CG14125	0	4.9471279	0.563	0.197	0
CG15818	0	3.2587424	0.492	0.108	0
Amy-d	0	3.4545206	0.718	0.233	0
CG15043	0	3.0710478	0.642	0.159	0
Mco1	0	2.4300353	0.424	0.061	0
CG13324	0	2.1388592	0.355	0.133	0
CG7589	0	2.2166688	0.386	0.075	0
Ag5r	0	2.4361222	0.307	0.066	0
CG13323	0	2.4128346	0.612	0.282	0
Adhr	0	1.6133663	0.541	0.202	0
CG42486	0	1.0222549	0.310	0.015	0
CG14933	0	1.4466203	0.368	0.045	0
CG18223	0	0.6894215	0.231	0.006	0
Adh	0	1.5147031	0.503	0.189	0
CG10405	0	1.0750603	0.246	0.027	0
lncRNA:CR40469	0	0.6702035	0.995	0.911	0
CG14439	0	0.7704120	0.241	0.010	0
CG42729	0	0.2725671	0.165	0.001	0
isoQC	0	1.4823857	0.358	0.228	0
CG9701	0	0.5324651	0.208	0.007	0

Cluster 10

markers_MAST10 = FindMarkers(integrated, ident.1=10, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST10, file.path("results", "markers_MAST", "cluster10-markers_MAST.csv"))
knitr::kable(head(markers_MAST10, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
18SrRNA-Psi:CR45861	0	1.9929792	0.774	0.425	0
18SrRNA:CR41548	0	1.8092135	0.917	0.516	0
18SrRNA:CR45838	0	1.8088839	0.917	0.516	0
18SrRNA:CR45841	0	1.8323186	0.902	0.506	0
pre-rRNA:CR45846	0	1.8543423	0.967	0.592	0
18SrRNA-Psi:CR41602	0	1.9522780	0.714	0.415	0
pre-rRNA:CR45856	0	1.9342579	0.842	0.497	0
pre-rRNA:CR45847	0	1.6122745	0.994	0.737	0
pre-rRNA:CR45845	0	1.5562434	0.979	0.810	0
28SrRNA-Psi:CR45851	0	1.3256885	0.688	0.531	0
28SrRNA-Psi:CR41609	0	1.3964064	0.560	0.459	0
28SrRNA:CR45844	0	1.1993126	0.744	0.607	0
28SrRNA-Psi:CR40741	0	1.3577011	0.565	0.448	0
28SrRNA:CR45837	0	0.7395038	0.411	0.486	0
CG12374	0	0.3401385	0.943	0.678	0
28SrRNA-Psi:CR45855	0	0.4141037	0.098	0.243	0
14-3-3epsilon	0	-0.6633088	0.158	0.488	0
Jon99Cii	0	0.5646288	0.893	0.650	0
Hr4	0	0.5679335	0.321	0.432	0
rg	0	0.3604389	0.137	0.328	0
Bace	0	0.2581049	0.923	0.674	0
His2Av	0	-0.3896197	0.062	0.330	0
Jon99Ciii	0	0.2742590	0.938	0.731	0
CG30025	0	0.2913271	0.899	0.671	0
CG42322	0	-0.2620372	0.182	0.471	0

Cluster 11

markers_MAST11 = FindMarkers(integrated, ident.1=11, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST11, file.path("results", "markers_MAST", "cluster11-markers_MAST.csv"))
knitr::kable(head(markers_MAST11, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Acbp3	0	2.7989831	0.978	0.298	0
CG32473	0	2.6193899	0.739	0.143	0
Gs2	0	2.5571275	0.978	0.369	0
CG15254	0	2.5365965	0.796	0.138	0
zetaTry	0	2.2494411	0.873	0.180	0
CG31266	0	2.3784007	0.672	0.081	0
NAAT1	0	1.9975076	0.669	0.058	0
CG8774	0	2.0324823	0.768	0.148	0
CG13492	0	1.8982504	0.876	0.217	0
CG31343	0	1.8074765	0.984	0.310	0
CG4053	0	1.8113458	0.573	0.030	0
Pdk	0	1.7664866	0.640	0.155	0
CG31267	0	2.2780699	0.564	0.030	0
CG31265	0	1.8975320	0.580	0.052	0
Gdh	0	1.7188263	0.768	0.192	0
Pepck2	0	1.7380777	0.653	0.096	0
CG9673	0	1.7049865	0.895	0.234	0
CG31269	0	1.3435526	0.433	0.007	0
CG31091	0	1.9780214	0.503	0.031	0
CG17475	0	1.9069711	0.411	0.005	0
CG32379	0	1.2844718	0.427	0.008	0
CG5958	0	1.7785181	0.726	0.177	0
CBP	0	0.2538569	0.344	0.001	0
Agpat4	0	1.5228590	0.761	0.194	0
CG10116	0	1.7645785	0.812	0.219	0

Cluster 12

markers_MAST12 = FindMarkers(integrated, ident.1=12, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST12, file.path("results", "markers_MAST", "cluster12-markers_MAST.csv"))
knitr::kable(head(markers_MAST12, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CCHa1	0	2.954678	0.603	0.119	0
IA-2	0	2.647187	0.993	0.403	0
nrv3	0	2.247964	0.909	0.242	0
Npc2b	0	2.621671	0.593	0.121	0
CG42541	0	2.090760	0.547	0.137	0
7B2	0	2.096422	0.915	0.231	0
w	0	2.291362	0.759	0.206	0
CCHa2	0	4.295649	0.407	0.249	0
Galphao	0	1.910813	0.691	0.175	0
pros	0	2.015971	0.769	0.181	0
Pal2	0	2.079188	0.782	0.147	0
CG42613	0	1.576964	0.450	0.043	0
CG30183	0	1.964119	0.736	0.130	0
Phm	0	1.793405	0.866	0.234	0
heph	0	1.849458	0.664	0.104	0
CG17646	0	1.608020	0.453	0.094	0
CG7191	0	2.818244	0.332	0.074	0
Tk	0	3.468051	0.557	0.218	0
Mrp4	0	1.856973	0.528	0.050	0
AstC	0	3.204312	0.524	0.302	0
amon	0	1.943279	0.638	0.127	0
Nrx-1	0	1.644528	0.570	0.095	0
CG32547	0	1.679896	0.463	0.056	0
unc-13-4A	0	2.278254	0.632	0.172	0
Hsp22	0	1.719189	0.720	0.342	0

Cluster 13

markers_MAST13 = FindMarkers(integrated, ident.1=13, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST13, file.path("results", "markers_MAST", "cluster13-markers_MAST.csv"))
knitr::kable(head(markers_MAST13, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CG42825	0	2.2269327	0.761	0.116	0
CG14499	0	3.0172332	0.764	0.075	0
lectin-37Da	0	2.7082242	0.757	0.081	0
Mur29B	0	2.5662759	0.804	0.184	0
CG33926	0	2.1060157	0.950	0.228	0
CG10912	0	2.2182013	0.907	0.203	0
CG32368	0	2.1692952	0.854	0.168	0
asRNA:CR44192	0	1.3426450	0.635	0.150	0
Npc2e	0	2.1977119	0.691	0.154	0
CG10911	0	1.5862509	0.963	0.401	0
CG31086	0	1.5386260	0.804	0.311	0
CG8353	0	1.6389178	0.535	0.155	0
CG31323	0	1.7430702	0.837	0.335	0
CG11912	0	2.4572814	0.412	0.086	0
CG10943	0	1.3910626	0.498	0.151	0
CG18493	0	1.4654841	0.754	0.253	0
CG9568	0	1.6351462	0.688	0.224	0
CG15255	0	1.3555371	0.771	0.267	0
lectin-37Db	0	1.3335771	0.688	0.226	0
RpL3	0	-0.9937845	0.718	0.896	0
CG13482	0	1.4957212	0.625	0.217	0
Cht4	0	0.7473134	0.309	0.016	0
Nazo	0	1.1559679	0.518	0.128	0
CG8950	0	0.3780558	0.455	0.065	0
CG14500	0	0.8257770	0.518	0.090	0

Cluster 14

markers_MAST14 = FindMarkers(integrated, ident.1=14, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST14, file.path("results", "markers_MAST", "cluster14-markers_MAST.csv"))
knitr::kable(head(markers_MAST14, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
NPF	0	5.1479016	0.983	0.290	0
chrb	0	3.5003289	0.841	0.281	0
IA-2	0	3.2135054	1.000	0.403	0
7B2	0	3.0448420	0.924	0.231	0
Phm	0	2.6654247	0.890	0.234	0
Tk	0	2.4481634	0.857	0.210	0
svr	0	2.5674166	0.867	0.243	0
Ih	0	2.4344804	0.661	0.097	0
CG15312	0	1.9534549	0.591	0.071	0
CG30183	0	2.1819966	0.764	0.129	0
nrv3	0	2.2510264	0.827	0.245	0
cpo	0	2.4463486	0.761	0.275	0
Hsp23	0	2.8419301	0.781	0.318	0
esg	0	2.2221630	0.761	0.184	0
Pal2	0	2.1567229	0.748	0.149	0
Hsp22	0	1.8522887	0.777	0.341	0
Sh	0	1.6991290	0.492	0.042	0
Hk	0	1.8939359	0.495	0.026	0
Ldh	0	2.2306744	0.728	0.177	0
CG46385	0	2.2108067	0.831	0.325	0
unc-13-4A	0	1.7157612	0.757	0.169	0
AstA-R2	0	1.3280768	0.425	0.014	0
Hsp70Aa	0	1.7015774	0.724	0.352	0
Or49b	0	0.8860271	0.342	0.002	0
Hsp70Ab	0	1.6730323	0.714	0.350	0

Cluster 15

markers_MAST15 = FindMarkers(integrated, ident.1=15, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST15, file.path("results", "markers_MAST", "cluster15-markers_MAST.csv"))
knitr::kable(head(markers_MAST15, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Orcokinin	0	5.3271710	1.000	0.292	0
AstC	0	3.3331658	0.996	0.293	0
CG4587	0	2.9631960	0.743	0.008	0
CG14989	0	3.7994074	0.871	0.190	0
CG13135	0	0.5614222	0.680	0.002	0
svr	0	3.0093252	0.979	0.244	0
CG33639	0	0.2836585	0.676	0.002	0
tap	0	1.9801462	0.851	0.030	0
RYa-R	0	0.3823188	0.593	0.000	0
unc-13-4A	0	2.9465008	0.963	0.167	0
Pvf3	0	0.3481559	0.676	0.005	0
CG30340	0	0.8985034	0.614	0.002	0
blot	0	0.6021470	0.680	0.009	0
Rgk3	0	2.1701331	0.768	0.022	0
Gad1	0	0.2640746	0.573	0.001	0
CG30116	0	1.2765806	0.714	0.013	0
dac	0	1.2441220	0.710	0.013	0
CG34458	0	0.7437293	0.685	0.010	0
Slob	0	0.7153288	0.685	0.011	0
brp	0	1.3592246	0.739	0.021	0
CG15203	0	0.8878283	0.544	0.001	0
CG15894	0	0.3864234	0.676	0.011	0
Frq1	0	0.5909662	0.614	0.009	0
CG3955	0	0.3744050	0.660	0.011	0
CG4341	0	0.5783358	0.685	0.014	0

Cluster 16

markers_MAST16 = FindMarkers(integrated, ident.1=16, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST16, file.path("results", "markers_MAST", "cluster16-markers_MAST.csv"))
knitr::kable(head(markers_MAST16, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
EbpIII	0	4.834704	0.920	0.223	0
spidey	0	3.329530	0.734	0.216	0
CG31522	0	3.284502	0.688	0.226	0
to	0	3.587573	0.726	0.064	0
Adhr	0	3.095023	0.819	0.200	0
Adh	0	3.083127	0.776	0.187	0
CG1124	0	3.623984	0.620	0.110	0
CG30197	0	3.532219	0.684	0.094	0
emp	0	2.438345	0.523	0.089	0
wat	0	3.231272	0.532	0.046	0
CG10237	0	2.613266	0.536	0.048	0
ple	0	2.881654	0.553	0.030	0
Msr-110	0	2.666640	0.612	0.116	0
CG31523	0	2.969051	0.658	0.148	0
ADPS	0	2.858588	0.473	0.050	0
Acbp2	0	2.549511	0.903	0.364	0
vir-1	0	2.554232	0.494	0.033	0
CG3097	0	2.313336	0.485	0.065	0
ATPCL	0	2.223677	0.527	0.170	0
Men	0	2.158755	0.515	0.191	0
CG8306	0	2.509803	0.494	0.060	0
CG43134	0	4.113872	0.354	0.116	0
Ace	0	1.772902	0.430	0.135	0
Cpr67B	0	2.731262	0.435	0.042	0
CG4660	0	2.353659	0.376	0.042	0

Cluster 17

markers_MAST17 = FindMarkers(integrated, ident.1=17, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST17, file.path("results", "markers_MAST", "cluster17-markers_MAST.csv"))
knitr::kable(head(markers_MAST17, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Skp2	0	5.3000419	1.000	0.389	0
CG34324	0	5.2816359	0.987	0.277	0
CG14645	0	5.1663525	1.000	0.282	0
CG34220	0	5.0936717	0.987	0.420	0
CG3906	0	5.0871549	0.991	0.217	0
Muc68D	0	5.0305282	0.933	0.207	0
CG11672	0	3.8428054	0.871	0.037	0
CG4783	0	3.2040502	0.956	0.160	0
Pgant4	0	2.4506775	0.662	0.014	0
Idgf4	0	2.4551618	0.782	0.072	0
CG30026	0	1.7591050	0.511	0.003	0
CG43673	0	2.9267167	0.467	0.020	0
opm	0	1.6150495	0.462	0.147	0
CG9988	0	1.2403834	0.351	0.001	0
CG31077	0	2.0244394	0.298	0.002	0
sgl	0	1.4856293	0.387	0.187	0
Cys	0	1.5106883	0.507	0.172	0
Pgant8	0	0.8396921	0.271	0.002	0
CG13810	0	0.9988078	0.329	0.008	0
CG45061	0	1.3525408	0.342	0.026	0
Gfat1	0	1.3953112	0.360	0.097	0
NijC	0	1.1420618	0.271	0.018	0
Ctl2	0	1.3710576	0.387	0.104	0
Hexo1	0	1.4742258	0.489	0.142	0
CG30025	0	-3.1632394	0.116	0.690	0

Cluster 18

markers_MAST18 = FindMarkers(integrated, ident.1=18, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST18, file.path("results", "markers_MAST", "cluster18-markers_MAST.csv"))
knitr::kable(head(markers_MAST18, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Npc2f	0	3.6101816	0.988	0.150	0
thetaTry	0	3.0983381	0.988	0.149	0
CG4830	0	2.6994819	0.789	0.017	0
CG18404	0	4.0870292	0.845	0.161	0
CG4563	0	2.4116003	0.770	0.061	0
Peritrophin-15a	0	2.6028515	0.950	0.184	0
CG9682	0	2.5655588	0.727	0.027	0
yip7	0	2.3389780	1.000	0.539	0
Jon65Aiv	0	2.1546865	1.000	0.694	0
Jon65Aiii	0	2.0314714	1.000	0.695	0
CG4734	0	1.8411515	0.826	0.132	0
CG30031	0	1.3828164	1.000	0.535	0
gammaTry	0	1.3828164	1.000	0.535	0
Bace	0	1.8774625	0.994	0.678	0
deltaTry	0	1.3582514	1.000	0.532	0
CG34040	0	0.5422801	0.640	0.082	0
LysB	0	1.3202083	0.826	0.197	0
28SrRNA-Psi:CR45853	0	0.3463674	0.422	0.025	0
Pebp1	0	1.5259964	0.944	0.592	0
CG10477	0	1.4562131	0.665	0.131	0
Pgcl	0	0.4049199	0.553	0.069	0
lncRNA:CR34335	0	-2.0779373	0.652	0.821	0
LysD	0	1.1953521	0.776	0.228	0
CG30025	0	1.0353096	1.000	0.673	0
LManII	0	0.9659594	0.708	0.172	0

Cluster 19

markers_MAST19 = FindMarkers(integrated, ident.1=19, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST19, file.path("results", "markers_MAST", "cluster19-markers_MAST.csv"))
knitr::kable(head(markers_MAST19, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CG6409	0	4.803699	0.966	0.145	0
CG42656	0	4.436801	0.910	0.082	0
IM33	0	3.965089	0.831	0.081	0
CG1143	0	2.846481	0.584	0.052	0
Irk2	0	2.881903	0.584	0.057	0
Ggamma30A	0	2.425736	0.584	0.061	0
tws	0	2.662732	0.640	0.170	0
CG14949	0	1.958402	0.315	0.027	0
CG16704	0	2.138345	0.360	0.030	0
Nha2	0	2.023512	0.382	0.042	0
Gdh	0	2.581828	0.685	0.205	0
CG9993	0	1.838703	0.371	0.012	0
Vha13	0	1.570255	0.978	0.679	0
CG18473	0	1.749039	0.337	0.018	0
Msr-110	0	2.274386	0.596	0.123	0
Gdap1	0	1.572230	0.270	0.149	0
Nop17l	0	1.449145	0.326	0.154	0
GLS	0	1.632297	0.180	0.046	0
lncRNA:CR40469	0	1.381869	1.000	0.913	0
bru1	0	1.433266	0.247	0.014	0
CG14933	0	1.864400	0.461	0.054	0
CG17999	0	1.509553	0.303	0.011	0
sesB	0	1.222895	0.978	0.775	0
Vha16-1	0	1.169516	0.978	0.791	0
CG7365	0	1.059547	0.225	0.015	0

Cluster 20

markers_MAST20 = FindMarkers(integrated, ident.1=20, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST20, file.path("results", "markers_MAST", "cluster20-markers_MAST.csv"))
knitr::kable(head(markers_MAST20, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
CG13285	0	3.8504561	0.963	0.020	0
CG44013	0	2.9005776	0.778	0.020	0
Msr-110	0	3.1788920	0.963	0.123	0
Idgf6	0	2.6238014	0.778	0.043	0
CG43394	0	2.6414936	0.685	0.028	0
GstD11	0	2.0589475	0.537	0.001	0
CG10096	0	2.8104611	0.667	0.012	0
CG5065	0	2.2587751	0.741	0.095	0
Idgf4	0	2.3408393	0.815	0.083	0
ELOVL	0	1.5613344	0.370	0.087	0
bond	0	2.0227798	0.481	0.016	0
Skp2	0	2.3129974	0.593	0.401	0
CG9336	0	1.7430299	0.500	0.041	0
CG17839	0	1.6639478	0.426	0.030	0
CG34220	0	2.2185041	0.556	0.432	0
yellow-d	0	1.4443748	0.333	0.002	0
Cys	0	1.8452262	0.611	0.177	0
trn	0	1.5692802	0.407	0.090	0
Gasp	0	1.5855129	0.389	0.008	0
Muc68D	0	1.8044905	0.519	0.221	0
pyr	0	1.1379466	0.296	0.205	0
CG34324	0	1.9788417	0.500	0.291	0
Baldspot	0	1.9471257	0.537	0.252	0
Snp	0	0.9650331	0.167	0.131	0
CG14645	0	2.2891081	0.593	0.296	0

Cluster 21

markers_MAST21 = FindMarkers(integrated, ident.1=21, test.use="MAST") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  arrange(p_val_adj)
readr::write_csv(markers_MAST21, file.path("results", "markers_MAST", "cluster21-markers_MAST.csv"))
knitr::kable(head(markers_MAST21, 25))

gene	p_val	avg_logFC	pct.1	pct.2	p_val_adj
Mur18B	0	3.763457	0.583	0.080	0
Alp2	0	3.229087	0.556	0.029	0
CG6726	0	3.195073	0.611	0.109	0
Alp4	0	4.009171	0.639	0.055	0
CG3168	0	3.775266	0.528	0.073	0
Spat	0	1.457579	0.278	0.026	0
Oatp58Dc	0	1.726650	0.222	0.015	0
CG2680	0	1.974845	0.361	0.026	0
CG10513	0	2.711481	0.361	0.037	0
CG7882	0	2.203185	0.167	0.016	0
CG10226	0	1.227024	0.250	0.032	0
Smvt	0	2.016667	0.306	0.015	0
CG10553	0	1.840484	0.333	0.016	0
CG7084	0	2.335300	0.361	0.091	0
CG6733	0	1.778015	0.306	0.028	0
CG10560	0	1.147329	0.139	0.020	0
Irk3	0	1.842469	0.333	0.018	0
CG31373	0	1.888380	0.333	0.112	0
CG10514	0	2.743189	0.333	0.009	0
Tps1	0	1.900322	0.222	0.009	0
HDAC6	0	2.516269	0.500	0.161	0
UGP	0	2.228555	0.500	0.160	0
CG13309	0	2.391852	0.333	0.008	0
CG14292	0	3.042500	0.417	0.098	0
CG11426	0	1.755369	0.306	0.109	0

saveRDS(integrated, file.path("results","integrated.rds"))